//#include <memory>
//#include <iostream>
//#include <fstream>
//#include <iomanip>
//
//#include <mpi.h>
//
//#include "SolverHypre.h"
//#include "Parallel/ParallelMPI.h"
//
//using namespace std;
//
//
//SolverHypre::SolverHypre(math::SMatrix* A_, math::Vector* b_, math::Vector* x_) :
//    A(A_),
//    b(b_),
//    x(x_)
//{
//
//    wall_time0 = new WallTime("wallTime0");
//    wall_time1 = new WallTime("wallTime1");
//    wall_time2 = new WallTime("wallTime2");
//    wall_time3 = new WallTime("wallTime3");
//    wall_time4 = new WallTime("wallTime4");
//
//
//}
//
//
//SolverHypre::~SolverHypre()
//{
//    //log_info<<"destroying solver";
//    finish();
//    //wall_time0->print();
//    //wall_time1->print();
//}
//
//
//void SolverHypre::operator()()
//{
//    //wall_time1->restart();
//    solve();
//    //wall_time1->update();
//}
//
//
//
//void SolverHypre::init()
//{
//    ParallelMPI* pmpi = config->pmpi;
//
//    num_procs = pmpi->get_size();
//    myid = pmpi->get_rank();
//
//    HYPRE_Init();
//    n = b->size();
//    solver_id = 0;
//    vis = 0;
//    print_system = 0;
//
//    log_info<<"hypre solver begin";
//    /* Preliminaries: want at least one processor per row */
//    if (n < num_procs) n = sqrt(num_procs) + 1;
//
//    /* Each processor knows only of its own rows - the range is denoted by ilower
//    and upper.  Here we partition the rows. We account for the fact that
//    N may not divide evenly by the number of processors. */
//    local_size = n/num_procs;
//    extra = n - local_size*num_procs;
//
//    ilower = local_size*myid;
//    ilower += min(myid, extra);
//
//    iupper = local_size*(myid+1);
//    iupper += min(myid+1, extra);
//    iupper = iupper - 1;
//
//    /* How many rows do I have? */
//    local_size = iupper - ilower + 1; 
//
//
//    log_info<<"hypre solver begin creating matrix";
//    /* Create the matrix.
//    Note that this is a square matrix, so we indicate the row partition
//    size twice (since number of rows = number of cols) */
//    HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A_hypre);
//
//    /* Choose a parallel csr format storage (see the User's Manual) */
//    HYPRE_IJMatrixSetObjectType(A_hypre, HYPRE_PARCSR);
//
//    /* Initialize before setting coefficients */
//    HYPRE_IJMatrixInitialize(A_hypre);
//    /* Now go through my local rows and set the matrix entries.
//    Each row has at most 5 entries. For example, if n=3:
//
//    A = [M -I 0; -I M -I; 0 -I M]
//    M = [4 -1 0; -1 4 -1; 0 -1 4]
//
//    Note that here we are setting one row at a time, though
//    one could set all the rows together (see the User's Manual).
//    */
//
//    int nnz;
//    /* OK to use constant-length arrays for CPUs
//    double values[nRowNnzMax];
//    int cols[nRowNnzMax];
//    nRowNnzMax: should be larger than max number of any row in A
//    */
//    double *values = (double *) malloc(30*sizeof(double));
//    int *cols = (int *) malloc(30*sizeof(int));
//    int *tmp = (int *) malloc(2*sizeof(int));
//
//    log_info<<"hypre solver begin set matrix";
//
//    //row byb row
//    /*
//    for (i = ilower; i <= iupper; i++)
//    {
//        nnz = 0;
//        if(i%(iupper/100) == 0) log_info<<"--------- progress = " << i/(iupper/100) << "%";
//
//        for(int j = 0; j < b->size(); j++)
//        {
//            double value;
//            bool status = A->get_value(i, j, value);
//            if(status)
//            {
//                cols[nnz] = j;
//                values[nnz] = value;
//                nnz++;
//            }
//        }
//
//        //Set the values for row i
//        tmp[0] = nnz;
//        tmp[1] = i;
//        HYPRE_IJMatrixSetValues(A_hypre, 1, &tmp[0], &tmp[1], cols, values);
//    }
//    */
//
//    //one by one
//    for (int j = 0; j < n; j++)
//    {
//        nnz = 0;
//        if(j%(iupper/10) == 0) log_info<<"--------- progress = " << j/(iupper/10) << "%";
//
//        auto& column_sm = A->data[j];
//        for(auto iter = column_sm.begin(); iter != column_sm.end(); iter++)
//        {
//            int I = iter->first;
//            if(I >= ilower && I <= iupper)
//            {
//                int J = j;
//                values[0] = iter->second;
//                cols[0] = J;
//                tmp[0] = 1;
//                tmp[1] = I;
//                HYPRE_IJMatrixSetValues(A_hypre, 1, &tmp[0], &tmp[1], cols, values);
//            }
//
//
//        }
//        
//    }
//
//    log_info<<"hypre solver end set matrix";
//
//    free(values);
//    free(cols);
//    free(tmp);
//
//
//    /* Assemble after setting the coefficients */
//    HYPRE_IJMatrixAssemble(A_hypre);
//
//    /* Note: for the testing of small problems, one may wish to read
//    in a matrix in IJ format (for the format, see the output files
//    from the -print_system option).
//    In this case, one would use the following routine:
//    HYPRE_IJMatrixRead( <filename>, MPI_COMM_WORLD,
//            HYPRE_PARCSR, &A );
//    <filename>  = IJ.A.out to read in what has been printed out
//    by -print_system (processor numbers are omitted).
//    A call to HYPRE_IJMatrixRead is an *alternative* to the
//    following sequence of HYPRE_IJMatrix calls:
//    Create, SetObjectType, Initialize, SetValues, and Assemble
//    */
//
//
//    /* Get the parcsr matrix object to use */
//    HYPRE_IJMatrixGetObject(A_hypre, (void**) &parcsr_A_hypre);
//
//
//    //log_info<<"hypre solver begin creating vector";
//    /* Create the rhs and solution */
//    HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&b_hypre);
//    HYPRE_IJVectorSetObjectType(b_hypre, HYPRE_PARCSR);
//    HYPRE_IJVectorInitialize(b_hypre);
//
//    HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&x_hypre);
//    HYPRE_IJVectorSetObjectType(x_hypre, HYPRE_PARCSR);
//    HYPRE_IJVectorInitialize(x_hypre);
//
//    rhs_values =  (double*) calloc(local_size, sizeof(double));
//    x_values =  (double*) calloc(local_size, sizeof(double));
//    rows = (int*) calloc(local_size, sizeof(int));
//
//    /* AMG */
//    if (solver_id == 0)
//    {
//        int num_iterations;
//        double final_res_norm;
//
//        /* Create solver */
//        HYPRE_BoomerAMGCreate(&solver);
//
//        /* Set some parameters (See Reference Manual for more parameters) */
//        //HYPRE_BoomerAMGSetPrintLevel(solver, 3);  /* print solve info + parameters */
//        HYPRE_BoomerAMGSetPrintLevel(solver, 0);  /* print solve info + parameters */
//        HYPRE_BoomerAMGSetOldDefault(solver); /* Falgout coarsening with modified classical interpolaiton */
//        HYPRE_BoomerAMGSetRelaxType(solver, config->em_param.solverRelaxType);   // G-S/Jacobi hybrid relaxation  3, 4, 6, 7, 18, 11, 12 
//        HYPRE_BoomerAMGSetRelaxOrder(solver, config->em_param.solverRelaxOrder);   /* uses C/F relaxation */
//        HYPRE_BoomerAMGSetNumSweeps(solver, config->em_param.solverNumSweeps);   /* Sweeeps on each level */
//        HYPRE_BoomerAMGSetMaxLevels(solver, config->em_param.solverMaxLevels);  // maximum number of levels, default 20
//        HYPRE_BoomerAMGSetTol(solver, config->em_param.solverRelativeTolerance);      /* conv. tolerance */
//
//
//        wall_time0->restart();
//        /* Set the rhs values and the solution to zero */
//        for (i=0; i<local_size; i++)
//        {
//            rhs_values[i] = (*b)(ilower+i);
//            x_values[i] = (*x)(ilower+i);
//            rows[i] = ilower + i;
//        }
//
//        HYPRE_IJVectorSetValues(b_hypre, local_size, rows, rhs_values);
//        HYPRE_IJVectorSetValues(x_hypre, local_size, rows, x_values);
//
//        wall_time0->update();
//
//
//        wall_time1->restart();
//        HYPRE_IJVectorAssemble(b_hypre);
//        /*  As with the matrix, for testing purposes, one may wish to read in a rhs:
//        HYPRE_IJVectorRead( <filename>, MPI_COMM_WORLD,
//                        HYPRE_PARCSR, &b );
//        as an alternative to the
//        following sequence of HYPRE_IJVectors calls:
//        Create, SetObjectType, Initialize, SetValues, and Assemble
//        */
//        HYPRE_IJVectorGetObject(b_hypre, (void **) &par_b_hypre);
//
//        HYPRE_IJVectorAssemble(x_hypre);
//        HYPRE_IJVectorGetObject(x_hypre, (void **) &par_x_hypre);
//        wall_time1->update();
//
//        /* Now setup and solve! */
//        wall_time2->restart();
//        HYPRE_BoomerAMGSetup(solver, parcsr_A_hypre, par_b_hypre, par_x_hypre);
//        wall_time2->update();
//
//    }
//
//
//}
//
//
//void SolverHypre::solve()
//{
//
//    if(config->sim_time->is_time_io())
//    {
//        HYPRE_BoomerAMGSetPrintLevel(solver, 3);
//    }
//    else
//    {
//        HYPRE_BoomerAMGSetPrintLevel(solver, 0);
//    }
//
//
//    //====================================update boundary===================================================
//    /*
//    HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A_hypre);
//
//    HYPRE_IJMatrixSetObjectType(A_hypre, HYPRE_PARCSR);
//
//    HYPRE_IJMatrixInitialize(A_hypre);
//
//
//    int nnz;
//
//    double *values = (double *) malloc(30*sizeof(double));
//    int *cols = (int *) malloc(30*sizeof(int));
//    int *tmp = (int *) malloc(2*sizeof(int));
//
//    log_info<<"hypre solver begin set matrix";
//    for (i = ilower; i <= iupper; i++)
//    {
//        nnz = 0;
//
//        for(int j = 0; j < b->size(); j++)
//        {
//            double value;
//            bool status = A->get_value(i, j, value);
//            if(status)
//            {
//                cols[nnz] = j;
//                values[nnz] = value;
//                nnz++;
//            }
//        }
//
//        tmp[0] = nnz;
//        tmp[1] = i;
//        HYPRE_IJMatrixSetValues(A_hypre, 1, &tmp[0], &tmp[1], cols, values);
//    }
//    log_info<<"hypre solver end set matrix";
//
//    free(values);
//    free(cols);
//    free(tmp);
//
//
//    HYPRE_IJMatrixAssemble(A_hypre);
//
//    HYPRE_IJMatrixGetObject(A_hypre, (void**) &parcsr_A_hypre);
//    */
//    //================================================================================================
//
//
//    wall_time0->restart();
//    /* Set the rhs values and the solution to zero */
//    for (i=0; i<local_size; i++)
//    {
//        rhs_values[i] = (*b)(ilower+i);
//        x_values[i] = (*x)(ilower+i);
//        rows[i] = ilower + i;
//    }
//
//    HYPRE_IJVectorSetValues(b_hypre, local_size, rows, rhs_values);
//    HYPRE_IJVectorSetValues(x_hypre, local_size, rows, x_values);
//
//    wall_time0->update();
//
//
//    wall_time1->restart();
//    HYPRE_IJVectorAssemble(b_hypre);
//    /*  As with the matrix, for testing purposes, one may wish to read in a rhs:
//    HYPRE_IJVectorRead( <filename>, MPI_COMM_WORLD,
//                    HYPRE_PARCSR, &b );
//    as an alternative to the
//    following sequence of HYPRE_IJVectors calls:
//    Create, SetObjectType, Initialize, SetValues, and Assemble
//    */
//    HYPRE_IJVectorGetObject(b_hypre, (void **) &par_b_hypre);
//
//    HYPRE_IJVectorAssemble(x_hypre);
//    HYPRE_IJVectorGetObject(x_hypre, (void **) &par_x_hypre);
//    wall_time1->update();
//
//    /*  Print out the system  - files names will be IJ.out.A.XXXXX
//    and IJ.out.b.XXXXX, where XXXXX = processor id */
//    if (print_system)
//    {
//        HYPRE_IJMatrixPrint(A_hypre, "IJ.out.A");
//        HYPRE_IJVectorPrint(b_hypre, "IJ.out.b");
//    }
//
//
//    /* Choose a solver and solve the system */
//    //log_info<<"hypre solver begin solving";
//
//    /* AMG */
//    if (solver_id == 0)
//    {
//        int num_iterations;
//        double final_res_norm;
//
//        /* Now setup and solve! */
//        wall_time2->restart();
//        //HYPRE_BoomerAMGSetup(solver, parcsr_A_hypre, par_b_hypre, par_x_hypre);
//        wall_time2->update();
//
//        wall_time3->restart();
//        HYPRE_BoomerAMGSolve(solver, parcsr_A_hypre, par_b_hypre, par_x_hypre);
//        wall_time3->update();
//
//
//        //log_info<<"hypre solver begin getting solution";
//        //get the solution
//        wall_time4->restart();
//        HYPRE_IJVectorGetValues(x_hypre, local_size, rows, x_values);
//        
//        //solution x now has only local solution
//        for(int i = 0; i < x->size(); i++)
//        {
//            if(i >= ilower && i <= iupper)
//            {
//                (*x)(i) = x_values[i-ilower];
//            }
//            else
//            {
//                (*x)(i) = 0.0;
//            }
//        }
//        wall_time3->update();    
//
//        /* Run info - needed logging turned on */
//        HYPRE_BoomerAMGGetNumIterations(solver, &num_iterations);
//        HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, &final_res_norm);
//        if (myid == 0 && print_system)
//        {
//            printf("\n");
//            printf("Iterations = %d\n", num_iterations);
//            printf("Final Relative Residual Norm = %e\n", final_res_norm);
//            printf("\n");
//        }
//        wall_time4->update();
//
//    }
//
//    //wall_time0->print();
//    //wall_time1->print();
//    //wall_time2->print();
//    //wall_time3->print();
//    //wall_time4->print();
//
//}
//
//void SolverHypre::finish()
//{
//    wall_time0->print();
//    wall_time1->print();
//    wall_time2->print();
//    wall_time3->print();
//
//    free(rhs_values);
//    free(rows);
//    free(x_values);
//
//    /* Destroy solver */
//    HYPRE_BoomerAMGDestroy(solver);
//}
